home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / bin / hinge / hinge.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-21  |  14.6 KB  |  706 lines
  1. #include <stdio.h>
  2. #include <signal.h>
  3. #include <math.h>
  4. #include "geom.h"
  5. #include "polylistP.h"
  6. #include "3d.h"
  7. #include "ooglutil.h"
  8. #include "hinge.h"
  9. #include "hui.h"
  10.  
  11. /*
  12.   pl is the polyhedron we're dealing with.
  13.  
  14.   int pl->n_verts:        total number of vertices in the polyhedron
  15.   HPoint3 pl->vl[]:        array containing vertices
  16.   int pl->n_polys:        total number of faces in the polyhedron
  17.   Poly pl->p[]:            array containing the faces
  18.   int pl->p[i].n_vertices:    number of vertices in i-th face
  19.   HPoint3 pl->p[i].v[j]->pt:    j-th vertex of i-th face
  20. */
  21. PolyList *pl;
  22.  
  23.  
  24. Transform *TT[MAXTHINGS];
  25.  
  26. float epsilon = .005;
  27.  
  28. /* pipes to and from geomview */
  29. FILE *togv, *fromgv;
  30.  
  31. char hingedatastring[] =
  32. #include "hingedata.h"
  33. ;
  34.  
  35. /* the space we're in */
  36. int space = HYPERBOLIC;
  37.  
  38. /* whether a rotation axis is currently defined */
  39. int haveaxis = 0;
  40.  
  41. /* the angle to rotate through */
  42. float angle = 90;
  43.  
  44. /* increments by which to rotate */
  45. float hingeincr = 5.0;
  46.  
  47. /* whether to kill geomview when we exit or die */
  48. int killgv = 0;
  49.  
  50. Point currentedge[2];
  51. Point base = {0, 0, 0, 1};
  52. Point tip = {0, 0, 1, 1};
  53. Point axis = {0, 0, 1, 1};
  54. Transform BaseT;
  55.  
  56. int mainpos[2];
  57. int helppos[2];
  58. int infopos[2];
  59. int filepos[2];
  60.  
  61. char *file = NULL;
  62.  
  63. static void space_alignZ(Point *a, Point *b, Transform T);
  64.  
  65. #define NEXTARG ++argv; --argc
  66.  
  67. usage()
  68. {
  69.   fprintf(stderr, "hinge [OPTIONS] [FILE]\n\
  70.     OPTIONS:\n\
  71.     -space {euclidean|hyperbolic}\n\
  72.     -angle a    (set angle initially to a degrees)\n\
  73.     -incr a        (hinge through increments of a degrees)\n\
  74.     -killgv        (kill geomview upon exit)\n\
  75.     -mpos x y    (position of main panel)\n\
  76.     -hpos x y    (position of help panel)\n\
  77.     -ipos x y    (position of info panel)\n\
  78.     -fpos x y    (position of file panel)\n");
  79.   exit(-1);
  80. }
  81.  
  82. void die(int sig, int code, struct sigcontext *sc)
  83. {
  84.   QuitProc(NULL,0);
  85. }
  86.  
  87. main(argc, argv)
  88.      int argc;
  89.      char *argv[];
  90. {
  91.   signal(SIGHUP,    die);
  92.   signal(SIGINT,    die);
  93.   signal(SIGQUIT,    die);
  94.   signal(SIGILL,    die);
  95.   signal(SIGTRAP,    die);
  96.   signal(SIGABRT,    die);
  97.   signal(SIGBUS,    die);
  98.   signal(SIGSEGV,    die);
  99.   signal(SIGPIPE,    die);
  100.   signal(SIGTERM,    die);
  101.   
  102.   NEXTARG;
  103.   
  104.   while (argc) {
  105.     if (!strcmp(argv[0],"-killgv")) {
  106.       killgv = 1;
  107.     } else if (!strcmp(argv[0],"-space")) {
  108.       NEXTARG;
  109.       if (!argc) usage();
  110.       if (!strcmp(argv[0],"hyperbolic")) {
  111.     space = HYPERBOLIC;
  112.       } else if (!strcmp(argv[0],"spherical")) {
  113.     space = SPHERICAL;
  114.       } else {;
  115.     space = EUCLIDEAN;
  116.       }
  117.     } else if (!strcmp(argv[0], "-incr")) {
  118.       NEXTARG;
  119.       if (!argc) usage();
  120.       hingeincr = atof(argv[0]);
  121.     } else if (!strcmp(argv[0], "-angle")) {
  122.       NEXTARG;
  123.       if (!argc) usage();
  124.       angle = atof(argv[0]);
  125.     } else if (!strcmp(argv[0], "-mpos")) {
  126.       NEXTARG;
  127.       if (argc<2) usage();
  128.       mainplacement = FL_PLACE_POSITION;
  129.       mainpos[0] = atoi(argv[0]);
  130.       mainpos[1] = atoi(argv[1]);
  131.       NEXTARG;
  132.     } else if (!strcmp(argv[0], "-fpos")) {
  133.       NEXTARG;
  134.       if (argc<2) usage();
  135.       fileplacement = FL_PLACE_POSITION;
  136.       filepos[0] = atoi(argv[0]);
  137.       filepos[1] = atoi(argv[1]);
  138.       NEXTARG;
  139.     } else if (!strcmp(argv[0], "-ipos")) {
  140.       NEXTARG;
  141.       if (argc<2) usage();
  142.       infoplacement = FL_PLACE_POSITION;
  143.       infopos[0] = atoi(argv[0]);
  144.       infopos[1] = atoi(argv[1]);
  145.       NEXTARG;
  146.     } else if (!strcmp(argv[0], "-hpos")) {
  147.       NEXTARG;
  148.       if (argc<2) usage();
  149.       helpplacement = FL_PLACE_POSITION;
  150.       helppos[0] = atoi(argv[0]);
  151.       helppos[1] = atoi(argv[1]);
  152.       NEXTARG;
  153.     } else {
  154.       file = argv[0];
  155.       if (argc>1) {
  156.     fprintf(stderr, "hinge: ignoring arguments after %s\n", file);
  157.     break;
  158.       }
  159.     }
  160.     
  161.     NEXTARG;
  162.   }
  163.   
  164.   hui_init();
  165.   Initialize();
  166.   hui_main_loop();
  167. }
  168.  
  169. Initialize()
  170. {
  171.   int i;
  172.  
  173.   LangInit();
  174.  
  175.   for (i=0; i<MAXTHINGS; ++i) TT[i] = NULL;
  176.   togv = stdout;
  177.   fromgv = stdin;
  178.   TmCopy(TM_IDENTITY, BaseT);
  179.   NewInst(0.0);
  180.   fprintf(togv, "(progn\n\
  181.   (normalization allgeoms none)\n\
  182.   (lines-closer allcams 1)\n\
  183.   (merge-baseap appearance { +edge })\n\
  184. )\n");
  185.   HingeSpace(space);
  186.   
  187.  
  188.   if (file == NULL) {
  189.     switch (space) {
  190.     case HYPERBOLIC:
  191.       file = "HingeDodec";
  192.       break;
  193.     default:
  194.     case EUCLIDEAN:
  195.       file = "HingeCube";
  196.       break;
  197.     }
  198.   }
  199.     
  200.   if (file != NULL) {
  201.     HingeLoad(file);
  202.   }
  203.  
  204.  
  205.   {
  206.     Geom *g = NULL;
  207.     FILE *fp = fopen("hingedata", "r");
  208.     if (fp == NULL)
  209.       fp = fstropen(hingedatastring, sizeof(hingedatastring), "r");
  210.     if (fp == NULL) {
  211.       OOGLError(0,"can't find the file \"hingedata\"; it must be\n\
  212. in the current directory.  This will be fixed soon. [mbp]");
  213.       exit(-1);
  214.     }
  215.     if (fp != NULL) {
  216.       g = GeomFLoad(fp, "hinge data");
  217.       fclose(fp);
  218.     }
  219.     if (g != NULL) {
  220.       DefinePick(g);
  221.     } else {
  222.       OOGLError(1,"can't read hinge data\n");
  223.       exit(-1);
  224.     }
  225.   }
  226.  
  227.   fprintf(togv, "(interest (pick world))\n");
  228.   fflush(togv);
  229. }
  230.  
  231. int WhichFace(HPoint3 *p, Transform T, PolyList *pl)
  232. {
  233.   int i;
  234.  
  235.   for (i=0; i<pl->n_polys; ++i) {
  236.     if (CoPlanar(p, T, &(pl->p[i]))) return i;
  237.   }
  238.   return -1;
  239. }
  240.  
  241. /*-----------------------------------------------------------------------
  242.  * Function:    CoPlanar
  243.  * Description:    test whether a point is coplanar with the image
  244.  *          of a polygon
  245.  * Args:    *p: the point
  246.  *        T: transform to apply to poly before test
  247.  *        *poly: the polygon
  248.  * Returns:    1 if yes, 0 if no
  249.  * Author:    mbp
  250.  * Date:    Thu Oct  1 10:06:39 1992
  251.  * Notes:    Uses global variable "epsilon" for test.
  252.  *        The plane of the polygon is assumed to be
  253.  *          that determined by its first 3 vertices.  This
  254.  *          assumes that those vertices are not collinear.
  255.  *          If they are collinear, this test may falsely
  256.  *          return 1.
  257.  *
  258.  */
  259. int CoPlanar(HPoint3 *p, Transform T, Poly *poly)
  260. {
  261.   register HPoint3 v0;
  262.   register HPoint3 v1;
  263.   register HPoint3 v2;
  264.   float det;
  265.   int ans;
  266.  
  267.   HPt3Transform(T, &(poly->v[0]->pt), &v0);
  268.   HPt3Transform(T, &(poly->v[1]->pt), &v1);
  269.   HPt3Transform(T, &(poly->v[2]->pt), &v2);
  270.  
  271.  
  272.   /* the following is the determinant of the 4x4 matrix
  273.      whose 1st 3 rows are the homog coords of the polygon's
  274.      1st 3 vertices, and whose last row is the homog coords
  275.      of p */
  276.   det = 
  277.     p->x * (  v2.w*v1.y*v0.z
  278.         - v1.w*v2.y*v0.z
  279.         - v2.w*v0.y*v1.z
  280.         + v0.w*v2.y*v1.z
  281.         + v1.w*v0.y*v2.z
  282.         - v0.w*v1.y*v2.z ) +
  283.     p->y * (- v2.w*v1.x*v0.z
  284.         + v1.w*v2.x*v0.z
  285.         + v2.w*v0.x*v1.z
  286.         - v0.w*v2.x*v1.z
  287.         - v1.w*v0.x*v2.z
  288.         + v0.w*v1.x*v2.z ) +
  289.     p->z * (  v2.w*v1.x*v0.y
  290.         - v1.w*v2.x*v0.y
  291.         - v2.w*v0.x*v1.y
  292.         + v0.w*v2.x*v1.y
  293.         + v1.w*v0.x*v2.y
  294.         - v0.w*v1.x*v2.y ) +
  295.     p->w * (- v2.x*v1.y*v0.z
  296.         + v1.x*v2.y*v0.z
  297.         + v2.x*v0.y*v1.z
  298.         - v0.x*v2.y*v1.z
  299.         - v1.x*v0.y*v2.z
  300.         + v0.x*v1.y*v2.z );
  301.  
  302.   ans = fabs(det) <= epsilon;
  303.   return ans;
  304. }
  305.  
  306.  
  307. float HPt3EucDistance( HPoint3 *a, HPoint3 *b )
  308. {
  309.   Point3 aa, ab;
  310.   float dist;
  311.  
  312.   HPt3ToPt3(a, &aa);
  313.   HPt3ToPt3(b, &ab);
  314.   dist = Pt3Distance( &aa, &ab );
  315.   return dist;
  316. }
  317.  
  318. int pt4equal(HPoint3 *a, HPoint3 *b)
  319. {
  320.   float dist;
  321.   int ans;
  322.  
  323.   dist = HPt3EucDistance( a, b );
  324.   ans = dist <= epsilon;
  325.   return ans;
  326. }
  327.  
  328. #if 0 /* commented out until spherical mode can be added; don't need
  329.      this now */
  330. float space_distance(HPoint3 *a, HPoint3 *b)
  331. {
  332.   float dist;
  333.   Point3 aa, bb;
  334.  
  335.   switch (space) {
  336.   default:
  337.   case EUCLIDEAN:
  338.     dist = HPt3EucDistance( a, b );
  339.     return dist;
  340.     break;
  341.   case HYPERBOLIC:
  342.     dist = HPt3HypDistance( a, b );
  343.     return dist;
  344.     break;
  345.   }
  346. }
  347.  
  348. float HPt3HypDistance( HPoint3 *a, HPoint3 *b )
  349. {
  350.   float ab, aa, bb, dist, p, s;
  351.   extern double    acosh(double);
  352.  
  353.   ab = MinkDot(a,b);
  354.   aa = MinkDot(a,a);
  355.   bb = MinkDot(b,b);
  356.   p  = ab*ab / (aa * bb);
  357.   if (p < 1) return 0;
  358.   s = sqrt(p);
  359.   dist = 2 * acosh( (double)s );
  360.   return dist;
  361. }
  362.  
  363. float MinkDot( HPoint3 *a, HPoint3 *b )
  364. {
  365.   return a->x*b->x + a->y*b->y + a->z*b->z - a->w*b->w;
  366. }
  367. #endif
  368.  
  369.  
  370. /*-----------------------------------------------------------------------
  371.  * Function:    PolyContainsEdge
  372.  * Description:    test whether the image of a polygon contain an edge
  373.  * Args:    e[]: array of two points --- the endpoints of the edge
  374.  *        T: transform to apply to the polygon
  375.  *        *poly: the polygon
  376.  * Returns:    -1, 0, or 1 (see below)
  377.  * Author:    mbp
  378.  * Date:    Thu Oct  1 12:28:36 1992
  379.  * Notes:    returns 1 if the polygon contains the edge with the vertices
  380.  *          in the same order, -1 if opposite order, 0 if polygon doesn't
  381.  *          contain the edge at all.
  382.  */
  383. int PolyContainsEdge(Point e[], Transform T, Poly *poly)
  384. {
  385.   int i,ans;
  386.   HPoint3 v;
  387.  
  388.   ans = 0;
  389.   for (i=0; i<poly->n_vertices; ++i) {
  390.     HPt3Transform(T, &(poly->v[i]->pt), &v);
  391.     if (pt4equal(&e[0], &v)) {
  392.       HPt3Transform(T, &(poly->v[(i+1)%poly->n_vertices]->pt), &v);
  393.       if (pt4equal(&e[1], &v)) {
  394.     ans = 1;
  395.     goto done;
  396.       }
  397.       HPt3Transform(T, &(poly->v[(i-1+poly->n_vertices)%poly->n_vertices]->pt), &v);
  398.       if (pt4equal(&e[1], &v)) {
  399.     ans = -1;
  400.     goto done;
  401.       }
  402.     }
  403.   }
  404.  done:
  405.   return ans;
  406. }
  407.  
  408. /* for debugging, just to make sure we understand what we have... */
  409. void WritePolyListInfo(PolyList *pl)
  410. {
  411.   int i,j;
  412.   FILE *fp = fopen("hinge.out", "w");
  413.  
  414.   fprintf(fp, "%1d vertices\n\n", pl->n_verts);
  415.   for (i=0; i<pl->n_verts; ++i) {
  416.     fprintf(fp, "vert[%2d] = (%6f, %6f, %6f, %6f)\n",
  417.         i,
  418.         pl->vl[i].pt.x,
  419.         pl->vl[i].pt.y,
  420.         pl->vl[i].pt.z,
  421.         pl->vl[i].pt.w);
  422.   }
  423.   fprintf(fp, "\n");
  424.   fprintf(fp, "%1d faces\n\n", pl->n_polys);
  425.   for (i=0; i<pl->n_polys; ++i) {
  426.     fprintf(fp, "face %1d has %1d vertices:\n", i, pl->p[i].n_vertices);
  427.     for (j=0; j<pl->p[i].n_vertices; ++j) {
  428.       fprintf(fp, "\t(%6f, %6f, %6f, %6f)\n",
  429.           pl->p[i].v[j]->pt.x,
  430.           pl->p[i].v[j]->pt.y,
  431.           pl->p[i].v[j]->pt.z,
  432.           pl->p[i].v[j]->pt.w);
  433.     }
  434.   }
  435.   fclose(fp);
  436. }
  437.  
  438. int
  439. HingeLoad(char *file)
  440. {
  441.   FILE *f;
  442.   Geom *g;
  443.  
  444.   if (strlen(file)<=0) return 0;
  445.  
  446.   /*
  447.     The following kludge uses geomview to get the data; we tell
  448.     geomview to load it, then dump it back to us, then delete it.
  449.    */
  450.   fprintf(togv,"(progn\n\
  451.   (geometry thing { < %s })\n\
  452.   (echo \"{\")\n\
  453.   (write geometry - thing self)\n\
  454.   (echo \"}\")\n\
  455.   (delete thing)\n\
  456. )\n", file);
  457.   fflush(togv);
  458.  
  459.   /*
  460.    * We now read the object we just told geomview to give us.
  461.    */
  462.   g = GeomFLoad(fromgv, "Hinge: pipe from geomview");
  463.  
  464.   if (g == NULL) {
  465.     fprintf(stderr, "Hinge: error reading geometry from file %s\n", file);
  466.   } else {
  467.     char buf[80];
  468.     DefineThing(g);
  469.     sprintf(buf,"%s",file);
  470.     hui_message(buf);
  471.     haveaxis = 0;
  472.     ShowAxis();
  473.   }
  474.   pl = (PolyList *)g;
  475.   return g != NULL;
  476. }
  477.   
  478.  
  479. int
  480. NewInst(float ang)
  481. {
  482.   int n = NewTTindex();
  483.   fprintf(togv,
  484.       "(progn (read geometry { define T%1d { LIST } } )\n\
  485. (geometry \"geom%1d\" { INST tlist { LIST { :T%1d } } geom :thing } ) )\n",
  486.       n, n, n);
  487.   fflush(togv);
  488.   Inst(n, 0.0);
  489.   return n;
  490. }
  491.  
  492. void
  493. ComputeTransform(Transform T, Transform BaseT, float ang)
  494. {
  495.   if (ang == 0) {
  496.     TmCopy(BaseT, T);
  497.   } else {
  498.     Rotation(T, ¤tedge[0], ¤tedge[1], RADIANS(ang));
  499.     TmConcat(BaseT, T, T);
  500.   }
  501. }
  502.  
  503. int
  504. NewTTindex()
  505. {
  506.   int i;
  507.   for (i=0; i<MAXTHINGS && TT[i]!=NULL; ++i);
  508.   if (i==MAXTHINGS) return -1;
  509. /*  TT[i] = (TmCoord *(*)[4])malloc(sizeof(Transform)); */
  510.   TT[i] = (void*)malloc(sizeof(Transform));
  511.   TmCopy(TM_IDENTITY, TT[i]);
  512.   return i;
  513. }
  514.  
  515. void
  516. DefineAxis(Point *a, Point *b)
  517. {
  518.   haveaxis = 1;
  519.   currentedge[0] = *a;
  520.   currentedge[1] = *b;
  521.   base = *a;
  522.   tip = *b;
  523.   HPt3Sub(b, a, &axis);
  524. }
  525.  
  526. void
  527. Rotation(Transform T, Point *a, Point *b, float angle)
  528. {
  529.   Transform A;
  530.  
  531.   space_alignZ(a, b, A);
  532.   TmInvert(A,A);
  533.   TmRotateZ(T, angle);
  534.   TmConjugate(T, A, T);
  535. }
  536.  
  537. void
  538. DefineThing(Geom *g)
  539. {
  540.   fprintf(togv,"(read geometry { define thing {\n");
  541.   GeomFSave(g, togv, "hinge output pipe (define thing...)");
  542.   fprintf(togv, "} } )\n");
  543.   fflush(togv);
  544. }
  545.  
  546. void
  547. DefinePick(Geom *g)
  548. {
  549.  
  550.   fprintf(togv, "(geometry \"pick\" { INST tlist { LIST { :edgeT } } geom { \n");
  551.   GeomFSave(g, togv, "hinge output pipe (geometry pick...)");
  552.   fprintf(togv, " } } )\n");
  553.  
  554.   fflush(togv);
  555. }
  556.  
  557. void
  558. Inst(int n, float ang)
  559. {
  560.   ComputeTransform(*(TT[n]), BaseT, ang);
  561.   fprintf(togv,"(read geometry { define T%1d TLIST\n",n);
  562.   fputtransform(togv, 1, (float*)(TT[n]), 0);
  563.   fprintf(togv," } )\n");
  564.   fflush(togv);
  565. }
  566.  
  567. void
  568. Undo()
  569. {
  570.   int i;
  571.  
  572.   for (i=0; i<MAXTHINGS && TT[i]!=NULL; ++i);
  573.   if (i==MAXTHINGS || i<=1) return;
  574.  
  575.   --i;
  576. /*  OOGLFree(TT[i]); */
  577.   free(TT[i]);
  578.   TT[i] = NULL;
  579.   fprintf(togv,"(delete \"geom%1d\")\n", i);
  580.   fflush(togv);
  581. }
  582.  
  583. void
  584. Reset()
  585. {
  586.   int i;
  587.  
  588.   fprintf(togv,"(progn\n");
  589.   for (i=1; i<MAXTHINGS && TT[i]!=NULL; ++i) {
  590.     /*  OOGLFree(TT[i]); */
  591.     free(TT[i]);
  592.     TT[i] = NULL;
  593.     fprintf(togv,"(delete \"geom%1d\")\n", i);
  594.   }
  595.   haveaxis = 0;
  596.   ShowAxis();
  597.   fprintf(togv,")\n");
  598.   fflush(togv);
  599. }
  600.  
  601. void
  602.   ShowAxis()
  603. {
  604.   Transform Ta, TaInv, R, Tloc, Tsize, Tnet;
  605.   Point3 a, b, b1;
  606.   float b1halflen, radius;
  607.   
  608.   if (haveaxis) {
  609.     
  610.     space_alignZ(¤tedge[0], ¤tedge[1], Tloc);
  611.     HPt3TransPt3(Tloc, ¤tedge[1], &b1);
  612.     TmInvert(Tloc, Tloc);
  613.     b1halflen = .5 * Pt3Length(&b1);
  614.     
  615.     TmTranslate( Tsize, (float)0, (float)0, b1halflen );
  616.     radius = .05 * b1halflen;
  617.     CtmScale( Tsize, radius, radius, b1halflen );
  618.     
  619.     TmConcat(Tsize, Tloc, Tnet);
  620.     
  621.     fprintf(togv,"(read geometry { define edgeT TLIST\n");
  622.     fputtransform(togv, 1, (float*)Tnet, 0);
  623.     fprintf(togv," } )\n");
  624.     
  625.   } else {
  626.     
  627.     fprintf(togv,"(read geometry { define edgeT { LIST } })\n");
  628.     
  629.   }
  630.   
  631.   fflush(togv);
  632. }
  633.  
  634. void space_translate_origin(Point *pt, Transform T)
  635. {
  636.   switch(space) {
  637.   case HYPERBOLIC:
  638.     TmHypTranslateOrigin(T, pt);
  639.     break;
  640.   case SPHERICAL:
  641.     TmSphTranslateOrigin(T, pt);
  642.     break;
  643.   default:
  644.   case EUCLIDEAN: TmTranslateOrigin(T, pt); break;
  645.   }
  646. }
  647.   
  648. void HingeSpace(int s)
  649. {
  650.   switch (s) {
  651.   case HYPERBOLIC:
  652.     space = HYPERBOLIC;
  653.     fprintf(togv, "\
  654. (progn\n\
  655.   (space hyperbolic)\n\
  656.   (camera-reset allcams)\n\
  657.   (merge-baseap appearance { +edge shading constant })\n\
  658. )\n");
  659.     fflush(togv);
  660.     break;
  661.   case SPHERICAL:
  662.     space = SPHERICAL;
  663.     fprintf(togv, "\
  664. (progn\n\
  665.   (space spherical)\n\
  666.   (camera-reset allcams)\n\
  667.   (merge-baseap appearance { +edge shading constant })\n\
  668. )\n");
  669.     fflush(togv);
  670.     break;
  671.   case EUCLIDEAN:
  672.   default:
  673.     space = EUCLIDEAN;
  674.     fprintf(togv, "\
  675. (progn\n\
  676.   (space euclidean)\n\
  677.   (merge-baseap appearance { +edge shading flat })\n\
  678. )\n");
  679.     fflush(togv);
  680.     break;
  681.   }
  682. }
  683.  
  684. /*
  685.   returns T = transform taking a to 0 and b to
  686.   a point on the +Z axis
  687.   */
  688. static void space_alignZ(Point *a, Point *b, Transform T)
  689. {
  690.   Transform Ta, R;
  691.   Point b1;
  692.   
  693.   /* Ta = transform taking a to origin */
  694.   space_translate_origin(a, Ta);
  695.   TmInvert(Ta,Ta);
  696.   
  697.   /* b1 = image of b under Ta */
  698.   HPt3Transform(Ta, b, &b1);
  699.   
  700.   /* R = transform rotating b1 to +Z axis */
  701.   TmRotateTowardZ(R, &b1);
  702.   
  703.   /* answer T is Ta * R */
  704.   TmConcat(Ta, R, T);
  705. }
  706.